## Purpose: analyzing depolicing data Philly

suppressPackageStartupMessages(
  {
    library(tidyverse)
    library(dplyr)
    library(estimatr)
    library(ggplot2)
    library(ggthemes)
    library(texreg)
    library(lfe)
    library(rdrobust)
    library(grid)
    library(gridExtra)
    library(rjson)
    library(tigris)
    library(sf)
    library(tidycensus)
  })

load("stopsveh_ph_geo.RData")
load("stopsped_ph_geo.RData")
load("crime_ph_geo.RData")

#### Creating rdit ####
## Crime
crime_ph_ts = 
  crime_ph_geo_ts %>% 
  mutate(bg = "fips") %>% 
  mutate(count = 1) %>% 
  group_by(bg, date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_ph_ts$blmrv = crime_ph_ts$date - min(crime_ph_ts$date[crime_ph_ts$blm == 1])

rdout_crime_ph = with(crime_ph_ts, rdrobust(x = blmrv, y = count, p = 1, all = TRUE))

## stops ped
stopsped_ph_ts = 
  stopsped_ph_geo_ts %>% 
  mutate(bg = "fips") %>% 
  mutate(count = 1) %>% 
  group_by(bg, date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsped_ph_ts$blmrv = stopsped_ph_ts$date - min(stopsped_ph_ts$date[stopsped_ph_ts$blm == 1])

rdout_stopsped_ph = with(stopsped_ph_ts, rdrobust(x = blmrv, y = count, p = 1, all = TRUE))

## stops veh
stopsveh_ph_ts = 
  stopsveh_ph_geo_ts %>% 
  mutate(bg = "fips") %>% 
  mutate(count = 1) %>% 
  group_by(bg, date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsveh_ph_ts$blmrv = stopsveh_ph_ts$date - min(stopsveh_ph_ts$date[stopsveh_ph_ts$blm == 1])

rdout_stopsveh_ph = with(stopsveh_ph_ts, rdrobust(x = blmrv, y = count, p = 1, all = TRUE))

##### Separate terciles to calculate total pop within these #####
## crime
crime_ph_geo3 = 
  crime_ph_geo_ts %>%
  mutate(inc_3 = ifelse(inc_low3 == 1, "low", "high")) %>%
  mutate(nw_3 = ifelse(nw_low3 == 1, "low", "high")) %>% 
  mutate(bl_3 = ifelse(bl_low3 == 1, "low", "high")) %>% 
  mutate(bg = as.numeric(fips))

crime_ph_geo_inc3 = 
  crime_ph_geo3 %>%
  group_by(bg) %>% 
  summarise(pop_total = mean(pop_total))

crime_ph_geo_inc3 = 
  crime_ph_geo3 %>%
  group_by(inc_3) %>% 
  summarise(pop_total_inc3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

crime_ph_geo_nw3 = 
  crime_ph_geo3 %>%
  group_by(nw_3) %>%
  summarise(pop_total_nw3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

crime_ph_geo_bl3 = 
  crime_ph_geo3 %>%
  group_by(bl_3) %>%
  summarise(pop_total_bl3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

crime_ph_geo3_ts = crime_ph_geo3
crime_ph_geo_inc3 <- as.data.frame(crime_ph_geo_inc3)
crime_ph_geo_nw3 <- as.data.frame(crime_ph_geo_nw3)
crime_ph_geo_bl3 <- as.data.frame(crime_ph_geo_bl3)
crime_ph_geo3_ts = merge(crime_ph_geo3_ts, crime_ph_geo_inc3, by = "inc_3")
crime_ph_geo3_ts = merge(crime_ph_geo3_ts, crime_ph_geo_nw3, by = "nw_3")
crime_ph_geo3_ts = merge(crime_ph_geo3_ts, crime_ph_geo_bl3, by = "bl_3")

## stops pedestrian
stopsped_ph_geo3 = 
  stopsped_ph_geo_ts %>%
  mutate(inc_3 = ifelse(inc_low3 == 1, "low", "high")) %>%
  mutate(nw_3 = ifelse(nw_low3 == 1, "low", "high")) %>% 
  mutate(bl_3 = ifelse(bl_low3 == 1, "low", "high")) %>% 
  mutate(bg = as.numeric(fips))

stopsped_ph_geo_inc3 = 
  stopsped_ph_geo3 %>%
  group_by(bg) %>% 
  summarise(pop_total = mean(pop_total))

stopsped_ph_geo_inc3 = 
  stopsped_ph_geo3 %>%
  group_by(inc_3) %>% 
  summarise(pop_total_inc3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

stopsped_ph_geo_nw3 = 
  stopsped_ph_geo3 %>%
  group_by(nw_3) %>%
  summarise(pop_total_nw3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

stopsped_ph_geo_bl3 = 
  stopsped_ph_geo3 %>%
  group_by(bl_3) %>%
  summarise(pop_total_bl3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

stopsped_ph_geo3_ts = stopsped_ph_geo3
stopsped_ph_geo_inc3 <- as.data.frame(stopsped_ph_geo_inc3)
stopsped_ph_geo_nw3 <- as.data.frame(stopsped_ph_geo_nw3)
stopsped_ph_geo_bl3 <- as.data.frame(stopsped_ph_geo_bl3)
stopsped_ph_geo3_ts = merge(stopsped_ph_geo3_ts, stopsped_ph_geo_inc3, by = "inc_3")
stopsped_ph_geo3_ts = merge(stopsped_ph_geo3_ts, stopsped_ph_geo_nw3, by = "nw_3")
stopsped_ph_geo3_ts = merge(stopsped_ph_geo3_ts, stopsped_ph_geo_bl3, by = "bl_3")

## stops vehicles
stopsveh_ph_geo3 = 
  stopsveh_ph_geo_ts %>%
  mutate(inc_3 = ifelse(inc_low3 == 1, "low", "high")) %>%
  mutate(nw_3 = ifelse(nw_low3 == 1, "low", "high")) %>% 
  mutate(bl_3 = ifelse(bl_low3 == 1, "low", "high")) %>% 
  mutate(bg = as.numeric(fips))

stopsveh_ph_geo_inc3 = 
  stopsveh_ph_geo3 %>%
  group_by(bg) %>% 
  summarise(pop_total = mean(pop_total))

stopsveh_ph_geo_inc3 = 
  stopsveh_ph_geo3 %>%
  group_by(inc_3) %>% 
  summarise(pop_total_inc3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

stopsveh_ph_geo_nw3 = 
  stopsveh_ph_geo3 %>%
  group_by(nw_3) %>%
  summarise(pop_total_nw3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

stopsveh_ph_geo_bl3 = 
  stopsveh_ph_geo3 %>%
  group_by(bl_3) %>%
  summarise(pop_total_bl3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

stopsveh_ph_geo3_ts = stopsveh_ph_geo3
stopsveh_ph_geo_inc3 <- as.data.frame(stopsveh_ph_geo_inc3)
stopsveh_ph_geo_nw3 <- as.data.frame(stopsveh_ph_geo_nw3)
stopsveh_ph_geo_bl3 <- as.data.frame(stopsveh_ph_geo_bl3)
stopsveh_ph_geo3_ts = merge(stopsveh_ph_geo3_ts, stopsveh_ph_geo_inc3, by = "inc_3")
stopsveh_ph_geo3_ts = merge(stopsveh_ph_geo3_ts, stopsveh_ph_geo_nw3, by = "nw_3")
stopsveh_ph_geo3_ts = merge(stopsveh_ph_geo3_ts, stopsveh_ph_geo_bl3, by = "bl_3")

#### Create plots ####
## crime
# income
crime_ph_geo_inclow3 = crime_ph_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "A. Crime (Income Lowest Tercile)")

crime_ph_geo_inchigh3 = crime_ph_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "B. Crime (Income Highest Tercile)")

# Nonwhite pop
crime_ph_geo_nwlow3 = crime_ph_geo3_ts %>% 
  filter(nw_low3 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "C. Crime (Nonwhite Population Lowest Tercile)")

crime_ph_geo_nwhigh3 = crime_ph_geo3_ts %>%
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "D. Crime (Nonwhite Population Highest Tercile)")

# Black pop
crime_ph_geo_bllow3 = crime_ph_geo3_ts %>% 
  filter(bl_low3 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "E. Crime (Black Population Lowest Tercile)")

crime_ph_geo_blhigh3 = crime_ph_geo3_ts %>%
  filter(bl_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "F. Crime (Black Population Highest Tercile)")

## Stops pedestrian
# income
stopsped_ph_geo_inclow3 = stopsped_ph_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "A. Pedestrian Stops (Income Lowest Tercile)")

stopsped_ph_geo_inchigh3 = stopsped_ph_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "B. Pedestrian Stops (Income Highest Tercile)")

# nonwhite pop
stopsped_ph_geo_nwlow3 = stopsped_ph_geo3_ts %>% 
  filter(nw_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "C. Pedestrian Stops (Nonwhite Pop Lowest Tercile)")

stopsped_ph_geo_nwhigh3 = stopsped_ph_geo3_ts %>% 
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "D. Pedestrian Stops (Nonwhite Pop Highest Tercile)")

# black pop
stopsped_ph_geo_bllow3 = stopsped_ph_geo3_ts %>% 
  filter(bl_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "E. Pedestrian Stops (Black Pop Lowest Tercile)")

stopsped_ph_geo_blhigh3 = stopsped_ph_geo3_ts %>% 
  filter(bl_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "F. Pedestrian Stops (Black Pop Highest Tercile)")

## Stops vehicles
# income
stopsveh_ph_geo_inclow3 = stopsveh_ph_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "A. Vehicle Stops (Income Lowest Tercile)")

stopsveh_ph_geo_inchigh3 = stopsveh_ph_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "B. Vehicle Stops (Income Highest Tercile)")

# nonwhite pop
stopsveh_ph_geo_nwlow3 = stopsveh_ph_geo3_ts %>% 
  filter(nw_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "C. Vehicle Stops (Nonwhite Pop Lowest Tercile)")

stopsveh_ph_geo_nwhigh3 = stopsveh_ph_geo3_ts %>% 
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "D. Vehicle Stops (Nonwhite Pop Highest Tercile)")

# black pop
stopsveh_ph_geo_bllow3 = stopsveh_ph_geo3_ts %>% 
  filter(bl_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "E. Vehicle Stops (Black Pop Lowest Tercile)")

stopsveh_ph_geo_blhigh3 = stopsveh_ph_geo3_ts %>% 
  filter(bl_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-30"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "F. Vehicle Stops (Black Pop Highest Tercile)")

## create and save grobs
crime_ph_tercileGrob = arrangeGrob(crime_ph_geo_inclow3, crime_ph_geo_inchigh3, 
                                crime_ph_geo_nwlow3, crime_ph_geo_nwhigh3, 
                                crime_ph_geo_bllow3, crime_ph_geo_blhigh3, 
                                ncol = 3, top=textGrob("Crime by Block Group (Philadelphia)"))

stopsped_ph_tercileGrob = arrangeGrob(stopsped_ph_geo_inclow3, stopsped_ph_geo_inchigh3,
                                stopsped_ph_geo_nwlow3, stopsped_ph_geo_nwhigh3,
                                stopsped_ph_geo_bllow3, stopsped_ph_geo_blhigh3,
                                ncol = 3, top=textGrob("Pedestrian Stops by Block Group (Philadelphia)"))

stopsveh_ph_tercileGrob = arrangeGrob(stopsveh_ph_geo_inclow3, stopsveh_ph_geo_inchigh3,
                                      stopsveh_ph_geo_nwlow3, stopsveh_ph_geo_nwhigh3,
                                      stopsveh_ph_geo_bllow3, stopsveh_ph_geo_blhigh3,
                                      ncol = 3, top=textGrob("Vehicle Stops by Block Group (Philadelphia)"))

ggsave(plot = crime_ph_tercileGrob, filename = "crime_geo_ph.png", width = 16, height = 8)
ggsave(plot = stopsped_ph_tercileGrob, filename = "stopsped_geo_ph.png", width = 16, height = 8)
ggsave(plot = stopsveh_ph_tercileGrob, filename = "stopsveh_geo_ph.png", width = 16, height = 8)

#### Analysis: RDs 25, 50, 100 bw
crime_ph_geo3_nwlow_25 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(nw_low3 == 1)

crime_ph_geo3_nwlow_50 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(nw_low3 == 1)

crime_ph_geo3_nwlow_100 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(nw_low3 == 1)

crime_ph_geo3_nwhigh_25 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(nw_high3 == 1)

crime_ph_geo3_nwhigh_50 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(nw_high3 == 1)

crime_ph_geo3_nwhigh_100 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(nw_high3 == 1)

crime_ph_geo3_inclow_25 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(inc_low3 == 1)

crime_ph_geo3_inclow_50 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(inc_low3 == 1)

crime_ph_geo3_inclow_100 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(inc_low3 == 1)

crime_ph_geo3_inchigh_25 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(inc_high3 == 1)

crime_ph_geo3_inchigh_50 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(inc_high3 == 1)

crime_ph_geo3_inchigh_100 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(inc_high3 == 1)

crime_ph_geo3_bllow_25 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(bl_low3 == 1)

crime_ph_geo3_bllow_50 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(bl_low3 == 1)

crime_ph_geo3_bllow_100 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(bl_low3 == 1)

crime_ph_geo3_blhigh_25 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(bl_high3 == 1)

crime_ph_geo3_blhigh_50 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(bl_high3 == 1)

crime_ph_geo3_blhigh_100 = crime_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(bl_high3 == 1)

stopsped_ph_geo3_nwlow_25 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(nw_low3 == 1)

stopsped_ph_geo3_nwlow_50 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(nw_low3 == 1)

stopsped_ph_geo3_nwlow_100 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(nw_low3 == 1)

stopsped_ph_geo3_nwhigh_25 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(nw_high3 == 1)

stopsped_ph_geo3_nwhigh_50 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(nw_high3 == 1)

stopsped_ph_geo3_nwhigh_100 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(nw_high3 == 1)

stopsped_ph_geo3_inclow_25 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(inc_low3 == 1)

stopsped_ph_geo3_inclow_50 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(inc_low3 == 1)

stopsped_ph_geo3_inclow_100 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(inc_low3 == 1)

stopsped_ph_geo3_inchigh_25 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(inc_high3 == 1)

stopsped_ph_geo3_inchigh_50 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(inc_high3 == 1)

stopsped_ph_geo3_inchigh_100 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(inc_high3 == 1)

stopsped_ph_geo3_bllow_25 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(bl_low3 == 1)

stopsped_ph_geo3_bllow_50 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(bl_low3 == 1)

stopsped_ph_geo3_bllow_100 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(bl_low3 == 1)

stopsped_ph_geo3_blhigh_25 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(bl_high3 == 1)

stopsped_ph_geo3_blhigh_50 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(bl_high3 == 1)

stopsped_ph_geo3_blhigh_100 = stopsped_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(bl_high3 == 1)

stopsveh_ph_geo3_nwlow_25 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(nw_low3 == 1)

stopsveh_ph_geo3_nwlow_50 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(nw_low3 == 1)

stopsveh_ph_geo3_nwlow_100 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(nw_low3 == 1)

stopsveh_ph_geo3_nwhigh_25 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(nw_high3 == 1)

stopsveh_ph_geo3_nwhigh_50 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(nw_high3 == 1)

stopsveh_ph_geo3_nwhigh_100 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(nw_high3 == 1)

stopsveh_ph_geo3_inclow_25 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(inc_low3 == 1)

stopsveh_ph_geo3_inclow_50 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(inc_low3 == 1)

stopsveh_ph_geo3_inclow_100 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(inc_low3 == 1)

stopsveh_ph_geo3_inchigh_25 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(inc_high3 == 1)

stopsveh_ph_geo3_inchigh_50 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(inc_high3 == 1)

stopsveh_ph_geo3_inchigh_100 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(inc_high3 == 1)

stopsveh_ph_geo3_bllow_25 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(bl_low3 == 1)

stopsveh_ph_geo3_bllow_50 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(bl_low3 == 1)

stopsveh_ph_geo3_bllow_100 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(bl_low3 == 1)

stopsveh_ph_geo3_blhigh_25 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-05-05") & date <= as.Date("2020-06-23")) %>% 
  filter(bl_high3 == 1)

stopsveh_ph_geo3_blhigh_50 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-04-09") & date <= as.Date("2020-07-18")) %>% 
  filter(bl_high3 == 1)

stopsveh_ph_geo3_blhigh_100 = stopsveh_ph_geo3_ts %>% 
  filter(date >= as.Date("2020-02-18") & date <= as.Date("2020-09-06")) %>% 
  filter(bl_high3 == 1)

## Separate by outcomes
crime_ph_geo3_inclow_ts = crime_ph_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_ph_geo3_inchigh_ts = crime_ph_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_ph_geo3_nwlow_ts = crime_ph_geo3_ts %>% 
  filter(nw_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_ph_geo3_nwhigh_ts = crime_ph_geo3_ts %>% 
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_ph_geo3_bllow_ts = crime_ph_geo3_ts %>% 
  filter(bl_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_ph_geo3_blhigh_ts = crime_ph_geo3_ts %>% 
  filter(bl_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsped_ph_geo3_inclow_ts = stopsped_ph_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsped_ph_geo3_inchigh_ts = stopsped_ph_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsped_ph_geo3_nwlow_ts = stopsped_ph_geo3_ts %>% 
  filter(nw_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/1000000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsped_ph_geo3_nwhigh_ts = stopsped_ph_geo3_ts %>% 
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsped_ph_geo3_bllow_ts = stopsped_ph_geo3_ts %>% 
  filter(bl_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsped_ph_geo3_blhigh_ts = stopsped_ph_geo3_ts %>% 
  filter(bl_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsveh_ph_geo3_inclow_ts = stopsveh_ph_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsveh_ph_geo3_inchigh_ts = stopsveh_ph_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsveh_ph_geo3_nwlow_ts = stopsveh_ph_geo3_ts %>% 
  filter(nw_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsveh_ph_geo3_nwhigh_ts = stopsveh_ph_geo3_ts %>% 
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsveh_ph_geo3_bllow_ts = stopsveh_ph_geo3_ts %>% 
  filter(bl_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stopsveh_ph_geo3_blhigh_ts = stopsveh_ph_geo3_ts %>% 
  filter(bl_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_ph_geo3_inclow_ts$blmrv = 
  crime_ph_geo3_inclow_ts$date - min(crime_ph_geo3_inclow_ts$date[crime_ph_geo3_inclow_ts$blm == 1])
crime_ph_geo3_inchigh_ts$blmrv = 
  crime_ph_geo3_inchigh_ts$date - min(crime_ph_geo3_inchigh_ts$date[crime_ph_geo3_inchigh_ts$blm == 1])
crime_ph_geo3_nwlow_ts$blmrv = 
  crime_ph_geo3_nwlow_ts$date - min(crime_ph_geo3_nwlow_ts$date[crime_ph_geo3_nwlow_ts$blm == 1])
crime_ph_geo3_nwhigh_ts$blmrv = 
  crime_ph_geo3_nwhigh_ts$date - min(crime_ph_geo3_nwhigh_ts$date[crime_ph_geo3_nwhigh_ts$blm == 1])
crime_ph_geo3_bllow_ts$blmrv = 
  crime_ph_geo3_bllow_ts$date - min(crime_ph_geo3_bllow_ts$date[crime_ph_geo3_bllow_ts$blm == 1])
crime_ph_geo3_blhigh_ts$blmrv = 
  crime_ph_geo3_blhigh_ts$date - min(crime_ph_geo3_blhigh_ts$date[crime_ph_geo3_blhigh_ts$blm == 1])

stopsped_ph_geo3_inclow_ts$blmrv = 
  stopsped_ph_geo3_inclow_ts$date - min(stopsped_ph_geo3_inclow_ts$date[stopsped_ph_geo3_inclow_ts$blm == 1])
stopsped_ph_geo3_inchigh_ts$blmrv = 
  stopsped_ph_geo3_inchigh_ts$date - min(stopsped_ph_geo3_inchigh_ts$date[stopsped_ph_geo3_inchigh_ts$blm == 1])
stopsped_ph_geo3_nwlow_ts$blmrv = 
  stopsped_ph_geo3_nwlow_ts$date - min(stopsped_ph_geo3_nwlow_ts$date[stopsped_ph_geo3_nwlow_ts$blm == 1])
stopsped_ph_geo3_nwhigh_ts$blmrv = 
  stopsped_ph_geo3_nwhigh_ts$date - min(stopsped_ph_geo3_nwhigh_ts$date[stopsped_ph_geo3_nwhigh_ts$blm == 1])
stopsped_ph_geo3_bllow_ts$blmrv = 
  stopsped_ph_geo3_bllow_ts$date - min(stopsped_ph_geo3_bllow_ts$date[stopsped_ph_geo3_bllow_ts$blm == 1])
stopsped_ph_geo3_blhigh_ts$blmrv = 
  stopsped_ph_geo3_blhigh_ts$date - min(stopsped_ph_geo3_blhigh_ts$date[stopsped_ph_geo3_blhigh_ts$blm == 1])

stopsveh_ph_geo3_inclow_ts$blmrv = 
  stopsveh_ph_geo3_inclow_ts$date - min(stopsveh_ph_geo3_inclow_ts$date[stopsveh_ph_geo3_inclow_ts$blm == 1])
stopsveh_ph_geo3_inchigh_ts$blmrv = 
  stopsveh_ph_geo3_inchigh_ts$date - min(stopsveh_ph_geo3_inchigh_ts$date[stopsveh_ph_geo3_inchigh_ts$blm == 1])
stopsveh_ph_geo3_nwlow_ts$blmrv = 
  stopsveh_ph_geo3_nwlow_ts$date - min(stopsveh_ph_geo3_nwlow_ts$date[stopsveh_ph_geo3_nwlow_ts$blm == 1])
stopsveh_ph_geo3_nwhigh_ts$blmrv = 
  stopsveh_ph_geo3_nwhigh_ts$date - min(stopsveh_ph_geo3_nwhigh_ts$date[stopsveh_ph_geo3_nwhigh_ts$blm == 1])
stopsveh_ph_geo3_bllow_ts$blmrv = 
  stopsveh_ph_geo3_bllow_ts$date - min(stopsveh_ph_geo3_bllow_ts$date[stopsveh_ph_geo3_bllow_ts$blm == 1])
stopsveh_ph_geo3_blhigh_ts$blmrv = 
  stopsveh_ph_geo3_blhigh_ts$date - min(stopsveh_ph_geo3_blhigh_ts$date[stopsveh_ph_geo3_blhigh_ts$blm == 1])

datasets_ph_geo = list(crime_ph_geo3_inclow_ts %>% as.data.frame %>% mutate(count_crimeph_inc_low = count) %>% filter(date >= as.Date("2016-01-01")),
                        crime_ph_geo3_inchigh_ts %>% as.data.frame %>% mutate(count_crimeph_inc_high = count) %>% filter(date >= as.Date("2016-01-01")),
                        crime_ph_geo3_nwlow_ts %>% as.data.frame %>% mutate(count_crimeph_nw_low = count) %>% filter(date >= as.Date("2016-01-01")),
                        crime_ph_geo3_nwhigh_ts %>% as.data.frame %>% mutate(count_crimeph_nw_high = count) %>% filter(date >= as.Date("2016-01-01")),
                        crime_ph_geo3_bllow_ts %>% as.data.frame %>% mutate(count_crimeph_bl_low = count) %>% filter(date >= as.Date("2016-01-01")),
                        crime_ph_geo3_blhigh_ts %>% as.data.frame %>% mutate(count_crimeph_bl_high = count) %>% filter(date >= as.Date("2016-01-01")),
                        stopsveh_ph_geo3_inclow_ts %>% as.data.frame %>% mutate(count_stopsvehph_inc_low = count) %>% filter(date >= as.Date("2016-01-01")),
                        stopsveh_ph_geo3_inchigh_ts %>% as.data.frame %>% mutate(count_stopsvehph_inc_high = count) %>% filter(date >= as.Date("2016-01-01")),
                        stopsveh_ph_geo3_nwlow_ts %>% as.data.frame %>% mutate(count_stopsvehph_nw_low = count) %>% filter(date >= as.Date("2016-01-01")),
                        stopsveh_ph_geo3_nwhigh_ts %>% as.data.frame %>% mutate(count_stopsvehph_nw_high = count) %>% filter(date >= as.Date("2016-01-01")),
                        stopsveh_ph_geo3_bllow_ts %>% as.data.frame %>% mutate(count_stopsvehph_bl_low = count) %>% filter(date >= as.Date("2016-01-01")),
                        stopsveh_ph_geo3_blhigh_ts %>% as.data.frame %>% mutate(count_stopsvehph_bl_high = count) %>% filter(date >= as.Date("2016-01-01")))

# outcomes

outcomes = c("count_crimeph_inc_low", "count_crimeph_inc_high", 
             "count_crimeph_nw_low", "count_crimeph_nw_high",
             "count_crimeph_bl_low", "count_crimeph_bl_high",
             "count_stopsvehph_inc_low", "count_stopsvehph_inc_high",
             "count_stopsvehph_nw_low", "count_stopsvehph_nw_high",
             "count_stopsvehph_bl_low", "count_stopsvehph_bl_high")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_ph_geo)))

for (i in 1:length(datasets_ph_geo)) {
  
  poly_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J Iteration ", j))
    
    kern_list = as.list(rep(NA, length(kerns)))
    
    for (l in 1:length(kerns)) {
      
      print(paste0("L Iteration ", l))
      
      mat_bw = matrix(NA, 100, 7)
      
      for(z in c(25, 50, 100)) {
        
        print(paste0("Z Iteration ", z))
        
        outrd = 
          rdrobust(x = as.numeric(datasets_ph_geo[[i]][, "blmrv"]),
                   y = datasets_ph_geo[[i]][, outcomes[i]],
                   p = polys[j],
                   h = z,
                   kernel = kerns[l])
        
        mat_bw[z, 1] = outrd$coef[3]
        mat_bw[z, 2] = outrd$se[3]
        mat_bw[z, 3] = outrd$pv[3]
        mat_bw[z, 4] = polys[j]
        mat_bw[z, 5] = kerns[l]
        mat_bw[z, 6] = outcomes[i]
        mat_bw[z, 7] = z
        
      }
      
      kern_list[[l]] = mat_bw %>%
        as.data.frame() %>%
        `colnames<-` (c("Coefficient", "SE", "Pvalue", "Polynomial", "Kernel", "Outcome", "BW"))%>% 
        mutate(Coefficient = as.numeric(as.character(Coefficient)),
               SE = as.numeric(as.character(SE)),
               Pvalue = as.numeric(as.character(Pvalue)),
               Polynomial = as.numeric(as.character(Polynomial)),
               Kernel = as.character(Kernel),
               Outcome = as.character(Outcome),
               BW = as.character(BW))
    }
    
    poly_list[[j]] = kern_list
    
  }
  
  dataset_list[[i]] = poly_list
  
}

ph_rdit_geo = dataset_list %>% 
  unlist(recursive = FALSE) %>%
  unlist(recursive = FALSE) %>%
  do.call(rbind.data.frame, .) %>% 
  mutate(Outcome = factor(Outcome,
                          levels = c("count_crimeph_inc_low", "count_crimeph_inc_high", 
                                     "count_crimeph_nw_low", "count_crimeph_nw_high",
                                     "count_crimeph_bl_low", "count_crimeph_bl_high",
                                     "count_stopspedph_inc_low", "count_stopspedph_inc_high",
                                     "count_stopspedph_nw_low", "count_stopspedph_nw_high",
                                     "count_stopspedph_bl_low", "count_stopspedph_bl_high",
                                     "count_stopsvehph_inc_low", "count_stopsvehph_inc_high",
                                     "count_stopsvehph_nw_low", "count_stopsvehph_nw_high",
                                     "count_stopsvehph_bl_low", "count_stopsvehph_bl_high"),
                          labels = c("A. Crime, Low Income",
                                     "B. Crime, High Income",
                                     "C. Crime, Low Percentage Nonwhite",
                                     "D. Crime, High Percentage Nonwhite",
                                     "E. Crime, Low Percentage Black",
                                     "F. Crime, High Percentage Black",
                                     "G. Pedestrian Stops, Low Income",
                                     "H. Pedestrian Stops, High Income",
                                     "I. Pedestrian Stops, Low Percentage Nonwhite",
                                     "J. Pedestrian Stops, High Percentage Nonwhite",
                                     "K. Pedestrian Stops, Low Percentage Black",
                                     "L. Pedestrian Stops, High Percentage Black",
                                     "G. Vehicle Stops, Low Income",
                                     "H. Vehicle Stops, High Income",
                                     "I. Vehicle Stops, Low Percentage Nonwhite",
                                     "J. Vehicle Stops, High Percentage Nonwhite",
                                     "K. Vehicle Stops, Low Percentage Black",
                                     "L. Vehicle Stops, High Percentage Black"
                          )),
         Polynomial = factor(Polynomial, levels = c(0, 1, 2),
                             labels = c("DIM", "Linear", "Quadratic")),
         Kernel = factor(Kernel, levels = c("tri", "uni", "epa"),
                         labels = c("Triangular",
                                    "Uniform",
                                    "Epanechnikov"))) %>% 
  filter(!is.na(BW))

ph_rdit_geo = ph_rdit_geo %>% 
  mutate(kern_poly = factor(paste0(Kernel, ", ", Polynomial), levels = c("Triangular, DIM", 
                                                                         "Triangular, Linear",
                                                                         "Triangular, Quadratic",
                                                                         "Uniform, DIM",
                                                                         "Uniform, Linear",
                                                                         "Uniform, Quadratic",
                                                                         "Epanechnikov, DIM",
                                                                         "Epanechnikov, Linear",
                                                                         "Epanechnikov, Quadratic")))

ph_rdit_geo$BW = as.numeric(ph_rdit_geo$BW)


#### Coefficient tests of difference ####
## Create function for s
s <- function(n1, n2, p, s1, s2){sqrt(((n1 - p)*((s1)^2) + (n2- p)*((s2)^2)/(n1 + n2 - 2*p)))}

se <- function(s, SEb1, SEb2, s1, s2){s*sqrt(((SEb1/s1)^2) + (SEb2/s2)^2)}

z <- function(b1, b2, se){(b1-b2)/se}

b_est <- function(b1, b2){(b1-b2)}

s_crimeph_25_inc <- s(n1 = 13169, n2 = 2648, p = 2, s1 = 0.0026084237, s2 = 0.0007399809)
se_cph_25_inc <- se(s = s_crimeph_25_inc, SEb1 = 0.0026084237, SEb2 = 0.0007399809, s1 = 0.0026084237, s2 = 0.0007399809)
z_cph_25_inc <- z(b1 = 2.455658e-03, b2 = 1.384564e-03, se = se_cph_25_inc)
b_cph_25_inc <- b_est(b1 = 2.455658e-03, b2 = 1.384564e-03)
p_cph_25_inc <- 2*pnorm(abs(z_cph_25_inc), mean = 0, lower.tail = FALSE)

s_crimeph_25_nw <- s(n1 = 267, n2 = 16878, p = 2, s1 = 0.0042557993, s2 = 0.0020423865)
se_cph_25_nw <- se(s = s_crimeph_25_nw, SEb1 = 0.0042557993, SEb2 = 0.0020423865, s1 = 0.0042557993, s2 = 0.0020423865)
z_cph_25_nw <- z(b1 = 5.534680e-03, b2 = 2.120341e-03, se = se_cph_25_nw)
b_cph_25_nw <- b_est(b1 = 5.534680e-03, b2 = 2.120341e-03)
p_cph_25_nw <- 2*pnorm(abs(z_cph_25_nw), mean = 0, lower.tail = FALSE)

s_crimeph_25_bl <- s(n1 = 839, n2 = 15605, p = 2, s1 = 0.0016924706, s2 = 0.0020592683)
se_cph_25_bl <- se(s = s_crimeph_25_bl, SEb1 = 0.0016924706, SEb2 = 0.0020592683, s1 = 0.0016924706, s2 = 0.0020592683)
z_cph_25_bl <- z(b1 = 3.755784e-03, b2 = 1.933889e-03, se = se_cph_25_bl)
b_cph_25_bl <- b_est(b1 = 3.755784e-03, b2 = 1.933889e-03)
p_cph_25_bl <- 2*pnorm(abs(z_cph_25_bl), mean = 0, lower.tail = FALSE)

s_crimeph_50_inc <- s(n1 = 23991, n2 = 5120, p = 2, s1 = 0.0015245411, s2 = 0.0004642689)
se_cph_50_inc <- se(s = s_crimeph_50_inc, SEb1 = 0.0015245411, SEb2 = 0.0004642689, s1 = 0.0015245411, s2 = 0.0004642689)
z_cph_50_inc <- z(b1 = 2.096633e-04, b2 = 7.575480e-04, se = se_cph_50_inc)
b_cph_50_inc <- b_est(b1 = 2.096633e-04, b2 = 7.575480e-04)
p_cph_50_inc <- 2*pnorm(abs(z_cph_50_inc), mean = 0, lower.tail = FALSE)

s_crimeph_50_nw <- s(n1 = 487, n2 = 30925, p = 2, s1 = 0.0027473960, s2 = 0.0011935714)
se_cph_50_nw <- se(s = s_crimeph_50_nw, SEb1 = 0.0027473960, SEb2 = 0.0011935714, s1 = 0.0027473960, s2 = 0.0011935714)
z_cph_50_nw <- z(b1 = 4.176682e-03, b2 = 3.804406e-04, se = se_cph_50_nw)
b_cph_50_nw <- b_est(b1 = 4.176682e-03, b2 = 3.804406e-04)
p_cph_50_nw <- 2*pnorm(abs(z_cph_50_nw), mean = 0, lower.tail = FALSE)

s_crimeph_50_bl <- s(n1 = 1531, n2 = 28705, p = 2, s1 = 0.0011170421, s2 = 0.0011962685)
se_cph_50_bl <- se(s = s_crimeph_50_bl, SEb1 = 0.0011170421, SEb2 = 0.0011962685, s1 = 0.0011170421, s2 = 0.0011962685)
z_cph_50_bl <- z(b1 = 2.754313e-03, b2 = 2.084095e-04, se = se_cph_50_bl)
b_cph_50_bl <- b_est(b1 = 2.754313e-03, b2 = 2.084095e-04)
p_cph_50_bl <- 2*pnorm(abs(z_cph_50_bl), mean = 0, lower.tail = FALSE)

s_crimeph_100_inc <- s(n1 = 48716, n2 = 10445, p = 2, s1 = 0.0008518542, s2 = 0.0002936343)
se_cph_100_inc <- se(s = s_crimeph_100_inc, SEb1 = 0.0008518542, SEb2 = 0.0002936343, s1 = 0.0008518542, s2 = 0.0002936343)
z_cph_100_inc <- z(b1 = -1.933809e-04, b2 = 2.572033e-04, se = se_cph_100_inc)
b_cph_100_inc <- b_est(b1 = -1.933809e-04, b2 = 2.572033e-04)
p_cph_100_inc <- 2*pnorm(abs(z_cph_100_inc), mean = 0, lower.tail = FALSE)

s_crimeph_100_nw <- s(n1 = 902, n2 = 62931, p = 2, s1 = 0.0017275033, s2 = 0.0006675944)
se_cph_100_nw <- se(s = s_crimeph_100_nw, SEb1 = 0.0017275033, SEb2 = 0.0006675944, s1 = 0.0017275033, s2 = 0.0006675944)
z_cph_100_nw <- z(b1 = 2.978004e-03, b2 = -5.803818e-05, se = se_cph_100_nw)
b_cph_100_nw <- b_est(b1 = 2.978004e-03, b2 = -5.803818e-05)
p_cph_100_nw <- 2*pnorm(abs(z_cph_100_nw), mean = 0, lower.tail = FALSE)

s_crimeph_100_bl <- s(n1 = 2972, n2 = 58453, p = 2, s1 = 0.0007676746, s2 = 0.0006649472)
se_cph_100_bl <- se(s = s_crimeph_100_bl, SEb1 = 0.0007676746, SEb2 = 0.0006649472, s1 = 0.0007676746, s2 = 0.0006649472)
z_cph_100_bl <- z(b1 = 2.148985e-03, b2 = -2.238355e-04, se = se_cph_100_bl)
b_cph_100_bl <- b_est(b1 = 2.148985e-03, b2 = -2.238355e-04)
p_cph_100_bl <- 2*pnorm(abs(z_cph_100_bl), mean = 0, lower.tail = FALSE)

s_stopsvehph_25_inc <- s(n1 = 4558, n2 = 364, p = 2, s1 = 0.0008980941, s2 = 0.0005273839)
se_svph_25_inc <- se(s = s_stopsvehph_25_inc, SEb1 = 0.0008980941, SEb2 = 0.0005273839, s1 = 0.0008980941, s2 = 0.0005273839)
z_svph_25_inc <- z(b1 = -9.553166e-03, b2 = -2.216620e-03, se = se_svph_25_inc)
b_svph_25_inc <- b_est(b1 = -9.553166e-03, b2 = -2.216620e-03)
p_svph_25_inc <- 2*pnorm(abs(z_svph_25_inc), mean = 0, lower.tail = FALSE)

s_stopsvehph_25_nw <- s(n1 = 23, n2 = 5388, p = 2, s1 = 0.0127822056, s2 = 0.0009507736)
se_svph_25_nw <- se(s = s_stopsvehph_25_nw, SEb1 = 0.0127822056, SEb2 = 0.0009507736, s1 = 0.0127822056, s2 = 0.0009507736)
z_svph_25_nw <- z(b1 = 1.465844e-02, b2 = -8.484643e-03, se = se_svph_25_nw)
b_svph_25_nw <- b_est(b1 = 1.465844e-02, b2 = -8.484643e-03)
p_svph_25_nw <- 2*pnorm(abs(z_svph_25_nw), mean = 0, lower.tail = FALSE)

s_stopsvehph_25_bl <- s(n1 = 84, n2 = 5238, p = 2, s1 = 0.0041109600, s2 = 0.0009063758)
se_svph_25_bl <- se(s = s_stopsvehph_25_bl, SEb1 = 0.0041109600, SEb2 = 0.0009063758, s1 = 0.0041109600, s2 = 0.0009063758)
z_svph_25_bl <- z(b1 = -2.718418e-04, b2 = -8.530001e-03, se = se_svph_25_bl)
b_svph_25_bl <- b_est(b1 = -2.718418e-04, b2 = -8.530001e-03)
p_svph_25_bl <- 2*pnorm(abs(z_svph_25_bl), mean = 0, lower.tail = FALSE)

s_stopsvehph_50_inc <- s(n1 = 8658, n2 = 634, p = 2, s1 = 0.0006248669, s2 = 0.0003625326)
se_svph_50_inc <- se(s = s_stopsvehph_50_inc, SEb1 = 0.0006248669, SEb2 = 0.0003625326, s1 = 0.0006248669, s2 = 0.0003625326)
z_svph_50_inc <- z(b1 = -9.679780e-03, b2 = -2.160428e-03, se = se_svph_50_inc)
b_svph_50_inc <- b_est(b1 = -9.679780e-03, b2 = -2.160428e-03)
p_svph_50_inc <- 2*pnorm(abs(z_svph_50_inc), mean = 0, lower.tail = FALSE)

s_stopsvehph_50_nw <- s(n1 = 36, n2 = 10261, p = 2, s1 = 0.0076236495, s2 = 0.0006037133)
se_svph_50_nw <- se(s = s_stopsvehph_50_nw, SEb1 = 0.0076236495, SEb2 = 0.0006037133, s1 = 0.0076236495, s2 = 0.0006037133)
z_svph_50_nw <- z(b1 = 4.874101e-03, b2 = -8.536182e-03, se = se_svph_50_nw)
b_svph_50_nw <- b_est(b1 = 4.874101e-03, b2 = -8.536182e-03)
p_svph_50_nw <- 2*pnorm(abs(z_svph_50_nw), mean = 0, lower.tail = FALSE)

s_stopsvehph_50_bl <- s(n1 = 169, n2 = 9923, p = 2, s1 = 0.0027555090, s2 = 0.0005847720)
se_svph_50_bl <- se(s = s_stopsvehph_50_bl, SEb1 = 0.0027555090, SEb2 = 0.0005847720, s1 = 0.0027555090, s2 = 0.0005847720)
z_svph_50_bl <- z(b1 = -8.130632e-04, b2 = -8.595830e-03, se = se_svph_50_bl)
b_svph_50_bl <- b_est(b1 = -8.130632e-04, b2 = -8.595830e-03)
p_svph_50_bl <- 2*pnorm(abs(z_svph_50_bl), mean = 0, lower.tail = FALSE)

s_stopsvehph_100_inc <- s(n1 = 27818, n2 = 2380, p = 2, s1 = 0.0005044425, s2 = 0.0002340851)
se_svph_100_inc <- se(s = s_stopsvehph_100_inc, SEb1 = 0.0005044425, SEb2 = 0.0002340851, s1 = 0.0005044425, s2 = 0.0002340851)
z_svph_100_inc <- z(b1 = -4.956212e-03, b2 = -6.994576e-04, se = se_svph_100_inc)
b_svph_100_inc <- b_est(b1 = -4.956212e-03, b2 = -6.994576e-04)
p_svph_100_inc <- 2*pnorm(abs(z_svph_100_inc), mean = 0, lower.tail = FALSE)

s_stopsvehph_100_nw <- s(n1 = 142, n2 = 33072, p = 2, s1 = 0.0042333885, s2 = 0.0004544779)
se_svph_100_nw <- se(s = s_stopsvehph_100_nw, SEb1 = 0.0042333885, SEb2 = 0.0004544779, s1 = 0.0042333885, s2 = 0.0004544779)
z_svph_100_nw <- z(b1 = 2.639920e-03, b2 = -4.385422e-03, se = se_svph_100_nw)
b_svph_100_nw <- b_est(b1 = 2.639920e-03, b2 = -4.385422e-03)
p_svph_100_nw <- 2*pnorm(abs(z_svph_100_nw), mean = 0, lower.tail = FALSE)

s_stopsvehph_100_bl <- s(n1 = 604, n2 = 31868, p = 2, s1 = 0.0018648278, s2 = 0.0004442453)
se_svph_100_bl <- se(s = s_stopsvehph_100_bl, SEb1 = 0.0018648278, SEb2 = 0.0004442453, s1 = 0.0018648278, s2 = 0.0004442453)
z_svph_100_bl <- z(b1 = 2.945763e-03, b2 = -4.501580e-03, se = se_svph_100_bl)
b_svph_100_bl <- b_est(b1 = 2.945763e-03, b2 = -4.501580e-03)
p_svph_100_bl <- 2*pnorm(abs(z_svph_100_bl), mean = 0, lower.tail = FALSE)
